TO APPEAR IN JCAP 



Use of event-level neutrino telescope 
data in global fits for theories of new 
physics 

P. Scott, 1,0 C. Savage, 2a J. Edsjo 2 and 
The IceCube Collaboration: 

R. Abbasi 3 , Y. Abdou 4 , M. Ackermann 5 , J. Adams 6 , J. A. Aguilar' 

M. Ahlers 3 , D. Altmann 8 , K. Andeen 3 , J. Auffenberg 3 , X. Bai 910 , 

M. Baker 3 , S. W. Barwick 11 , V. Baum 12 , R. Bay 13 , K. Beattie 14 , 

J. J. Beatty 15 16 , S. Bechet 17 , J. Becker Tjus 18 , K.-H. Becker 19 , 

M. Bell 20 , M. L. Benabderrahmane 5 , S. BenZvi 3 , J. Berdermann 5 , 

P. Berghaus 5 , D. Berley 21 , E. Bernardini 5 , D. Bertrand 17 , 

D. Z. Besson 22 , D. Bindig 19 , M. Bissok 23 , E. Blaufuss 21 , 

J. Blumenthal 23 , D. J. Boersma 23 , C. Bohm 2 , D. Bose 24 , 

S. Boser 25 , 0. Botner 26 , L. Brayeur 24 , A. M. Brown 6 , R. Bruijn 27 , 

J. Brunner 5 , S. Buitink 24 , K. S. Caballero-Mora 20 , M. Carson 4 , 

J. Casey 28 , M. Casier 24 , D. Chirkin 3 , B. Christy 21 , F. Clevermann 2 ^ 

S. Cohen 27 , D. F. Cowen 20 30 , A. H. Cruz Silva 5 , M. Danninger 2a , 

J. Daughhetee 28 , J. C. Davis 15 , C. De Clercq 24 , F. Descamps 3 , 

P. Desiati 3 , G. de Vries-Uiterweerd 4 , T. DeYoung 20 , 

J. C. Di'az-Velez 3 , J. Dreyer 18 , J. P. Dumm 3 , M. Dunkman 20 , 

R. Eagan 20 , J. Eisch 3 , R. W. Ellsworth 21 , O. Engdegard 26 , 

S. Euler 23 , P. A. Evenson 9 , O. Fadiran 3 , A. R. Fazely 31 , 

A. Fedynitch 18 , J. Feintzeig 3 , T. Feusels 4 , K. Filimonov 13 , 

C. Finley 2 , T. Fischer-Wasels 19 , S. Flis 2 , A. Franckowiak 25 , 

R. Franke 5 , K. Frantzen 29 , T. Fuchs 29 , T. K. Gaisser 9 , 

J. Gallagher 32 , L. Gerhardt 14 13 , L. Gladstone 3 , T. Glusenkamp 5 , 

A. Goldschmidt 14 , J. A. Goodman 21 , D. Gora 5 , D. Grant 33 , 

A. GroB 34 , S. Grullon 3 , M. Gurtner 19 , C. Ha 1413 , A. Haj Ismail 4 , 

A. Hallgren 26 , F. Halzen 3 , K. Hanson 17 , D. Heereman 17 , 

P. Heimann 23 , D. Heinen 23 , K. Helbing 19 , R. Hellauer 21 , 

S. Hickford 6 , G. C. Hill 35 , K. D. Hoffman 21 , R. Hoffmann 19 , 

A. Homeier 25 , K. Hoshina 3 , W. Huelsnitz 21 36 , P. O. Hulth 2 , 



K. Hultqvist 2 , S. Hussain 9 , A. Ishihara 37 , E. Jacobi 5 , J. Jacobsen 3 , 

G. S. Japaridze 38 , O. Jlelati 4 , H. Johansson 2 , A. Kappes 8 , 
T. Karg 5 , A. Karle 3 , J. Kiryluk 39 , F. Kislat 5 , J. Klas 19 , 

S. R. Klein 1413 , J.-H. Kohne 29 , G. Kohnen 40 , H. Kolanoski 8 , 
L. Kopke 12 , C. Kopper 3 , S. Kopper 19 , D. J. Koskinen 20 , 
M. Kowalski 25 , M. Krasberg 3 , G. Kroll 12 , J. Kunnen 24 , 
N. Kurahashi 3 , T. Kuwabara 9 , M. Labare 24 , K. Laihem 23 , 

H. Landsman 3 , M. J. Larson 41 , R. Lauer 5 , M. Lesiak-Bzdak 39 , 
J. Lunemann 12 , J. Madsen 42 , R. Maruyama 3 , K. Mase 37 , 

H. S. Matis 14 , F. McNally 3 , K. Meagher 21 , M. Merck 3 , 

P. Meszaros 30 20 , T. Meures 17 , S. Miarecki 14 13 , E. Middell 5 , 

N. Milke 29 , J. Miller 24 , L. Mohrmann 5 , T. Montaruli 743 , R. Morse 3 , 

S. M. Movit 30 , R. Nahnhauer 5 , U. Naumann 19 , S. C. Nowicki 33 , 

D. R. Nygren 14 , A. Obertacke 19 , S. Odrowski 34 , A. Olivas 21 , 

M. Olivo 18 , A. O'Murchadha 17 , S. Panknin 25 , L. Paul 23 , 

J. A. Pepper 41 , C. Perez de los Heros 26 , D. Pieloth 29 , N. Pirk 5 , 

J. Posselt 19 , P. B. Price 13 , G. T. Przybylski 14 , L. Radel 23 , 

K. Rawlins 44 , P. Redl 21 , E. Resconi 34 , W. Rhode 29 , M. Ribordy 27 , 

M. Richman 21 , B. Riedel 3 , J. P. Rodrigues 3 , F. Rothmaier 12 , 

C. Rott 15 , T. Ruhe 29 , D. Rutledge 20 , B. Ruzybayev 9 , 

D. Ryckbosch 4 , S. M. Saba 18 , T. Salameh 20 , H.-G. Sander 12 , 
M. Santander 3 , S. Sarkar 45 , K. Schatto 12 , M. Scheel 23 , 

F. Scheriau 29 , T. Schmidt 21 , M. Schmitz 29 , S. Schoenen 23 , 

S. Schoneberg 18 , L. Schonherr 23 , A. Schonwald 5 , A. Schukraft 23 , 
L. Schulte 25 , O. Schulz 34 , D. Seckel 9 , S. H. Seo 2 , Y. Sestayo 34 , 
S. Seunarine 46 , M. W. E. Smith 20 , M. Soiron 23 , D. Soldin 19 , 

G. M. Spiczak 42 , C. Spiering 5 , M. Stamatikos 15 47 , T. Stanev 9 , 
A. Stasik 25 , T. Stezelberger 14 , R. G. Stokstad 14 , A. StoBI 5 , 

E. A. Strahler 24 , R. Strom 26 , G. W. Sullivan 21 , H. Taavola 26 , 

I. Taboada 28 , A. Tamburro 9 , S. Ter-Antonyan 31 , S. Tilav 9 , 

P. A. Toale 41 , S. Toscano 3 , M. Usner 25 , N. van Eijndhoven 24 , 

D. van der Drift 1413 , A. Van Overloop 4 , J. van Santen 3 , 

M. Vehring 23 , M. Voge 25 , C. Walck 2 , T. Waldenmaier 8 , 

M. Wallraff 23 , M. Walter 5 , R. Wasserman 20 , Ch. Weaver 3 , 

C. Wendt 3 , S. Westerhoff 3 , N. Whitehorn 3 , K. Wiebe 12 , 

C. H. Wiebusch 23 , D. R. Williams 41 , H. Wissing 21 , M. Wolf 2 , 

T. R. Wood 33 , K. Woschnagg 13 , C. Xu 9 , D. L. Xu 41 , X. W. Xu 31 , 

J. P. Yanez 5 , G. Yodh 11 , S. Yoshida 37 , P. Zarzhitsky 41 , 

J. Ziemann 29 , A. Zilles 23 , and M. Zoll 2 

"Corresponding authors: P. Scott, M. Danninger, C. Savage 



x Dept. of Physics, McGill University, 3600 rue University, Montreal, QC, H3A 2T8, Canada 
2 Oskar Klein Centre for Cosmoparticle Physics and Dept. of Physics, Stockholm University, 
SE-10691 Stockholm, Sweden 

3 Dept. of Physics and Wisconsin IceCube Particle Astrophysics Center, University of Wis- 
consin, Madison, WI 53706, USA 
4 Dept. of Physics and Astronomy, University of Gent, B-9000 Gent, Belgium 
5 DESY, D-15735 Zeuthen, Germany 

6 Dept. of Physics and Astronomy, University of Canterbury, Private Bag 4800, Christchurch, 
New Zealand 

7 Departement de physique nucleaire et corpusculaire, Universite de Geneve, CH-1211 Geneve, 
Switzerland 

8 Institut fiir Physik, Humboldt-Universitat zu Berlin, D-12489 Berlin, Germany 

9 Bartol Research Institute and Department of Physics and Astronomy, University of Delaware, 

Newark, DE 19716, USA 
10 Physics Department, South Dakota School of Mines and Technology, Rapid City, SD 57701, 

USA 

11 Dept. of Physics and Astronomy, University of California, Irvine, CA 92697, USA 
12 Institute of Physics, University of Mainz, Staudinger Weg 7, D-55099 Mainz, Germany 
13 Dept. of Physics, University of California, Berkeley, CA 94720, USA 
14 Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA 

15 Dept. of Physics and Center for Cosmology and Astro-Particle Physics, Ohio State Univer- 
sity, Columbus, OH 43210, USA 
16 Dept. of Astronomy, Ohio State University, Columbus, OH 43210, USA 
17 Universite Libre de Bruxelles, Science Faculty CP230, B-1050 Brussels, Belgium 
18 Fakultat fiir Physik & Astronomie, Ruhr-Universitat Bochum, D-44780 Bochum, Germany 
19 Dept. of Physics, University of Wuppertal, D-42119 Wuppertal, Germany 
20 Dept. of Physics, Pennsylvania State University, University Park, PA 16802, USA 
21 Dept. of Physics, University of Maryland, College Park, MD 20742, USA 
22 Dept. of Physics and Astronomy, University of Kansas, Lawrence, KS 66045, USA 
23 III. Physikalisches Institut, RWTH Aachen University, D-52056 Aachen, Germany 
24 Vrije Universiteit Brussel, Dienst ELEM, B-1050 Brussels, Belgium 
25 Physikalisches Institut, Universitat Bonn, Nussallee 12, D-53115 Bonn, Germany 
26 Dept. of Physics and Astronomy, Uppsala University, Box 516, S-75120 Uppsala, Sweden 
27 Laboratory for High Energy Physics, Ecole Polytechnique Federale, CH-1015 Lausanne, 
Switzerland 

28 School of Physics and Center for Relativistic Astrophysics, Georgia Institute of Technology, 

Atlanta, GA 30332, USA 
29 Dept. of Physics, TU Dortmund University, D-44221 Dortmund, Germany 
30 Dept. of Astronomy and Astrophysics, Pennsylvania State University, University Park, PA 

16802, USA 

31 Dept. of Physics, Southern University, Baton Rouge, LA 70813, USA 
32 Dept. of Astronomy, University of Wisconsin, Madison, WI 53706, USA 
33 Dept. of Physics, University of Alberta, Edmonton, Alberta, Canada T6G 2G7 
34 T.U. Munich, D-85748 Garching, Germany 



35 School of Chemistry & Physics, University of Adelaide, Adelaide SA, 5005 Australia 
36 Los Alamos National Laboratory, Los Alamos, NM 87545, USA 
37 Dept. of Physics, Chiba University, Chiba 263-8522, Japan 
38 CTSPS, Clark-Atlanta University, Atlanta, GA 30314, USA 

39 Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794- 
3800, USA 

40 Universite de Mons, 7000 Mons, Belgium 

41 Dept. of Physics and Astronomy, University of Alabama, Tuscaloosa, AL 35487, USA 
42 Dept. of Physics, University of Wisconsin, River Falls, WI 54022, USA 
43 also Sezione INFN, Dipartimento di Fisica, 1-70126, Bari, Italy 

44 Dept. of Physics and Astronomy, University of Alaska Anchorage, 3211 Providence Dr., 

Anchorage, AK 99508, USA 
45 Dept. of Physics, University of Oxford, 1 Keble Road, Oxford OX1 3NP, UK 
46 Dept. of Physics, University of the West Indies, Cave Hill Campus, Bridgetown BB11000, 

Barbados 

47 NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA 
E-mail: patscott@physics.mcgill.ca, danning@fysik.su.se, savage@fysik.su.se 

Abstract. We present a fast likelihood method for including event-level neutrino telescope 
data in parameter explorations of theories for new physics, and announce its public release 
as part of DarkSUSY 5.0.6. Our construction includes both angular and spectral informa- 
tion about neutrino events, as well as their total number. We also present a corresponding 
measure for simple model exclusion, which can be used for single models without reference 
to the rest of a parameter space. We perform a number of supersymmetric parameter scans 
with IceCube data to illustrate the utility of the method: example global fits and a signal 
recovery in the constrained minimal supersymmetric standard model (CMSSM), and a model 
exclusion exercise in a 7-parameter phenomenological version of the MSSM. The final Ice- 
Cube detector configuration will probe almost the entire focus-point region of the CMSSM, 
as well as a number of MSSM-7 models that will not otherwise be accessible to e.g. direct 
detection. Our method accurately recovers the mock signal, and provides tight constraints on 
model parameters and derived quantities. We show that the inclusion of spectral information 
significantly improves the accuracy of the recovery, providing motivation for its use in future 
IceCube analyses. 

Keywords: dark matter theory, dark matter experiments, neutrino astronomy, cosmology 
of theories beyond the SM 
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1 Introduction 

Despite ongoing efforts, we have yet to identify dark matter. One of the most promising 
candidate classes is the weakly-interacting massive particle (WIMP) [1-3], so named because 
WIMPs interact with standard model (SM) particles only via the weak nuclear force. WIMPs 
are appealing because the expected cosmological density of a thermal relic with a weak-scale 
annihilation cross-section is the same order of magnitude as the observed value [4]. 

WIMPs are typically sought via three observational channels: direct WIMP-nucleon 
scattering (e.g. [5-11]), production at accelerators (e.g. [12-18]) and indirect detection of SM 
products of WIMP self-annihilation (e.g. [19-25]). A promising version of indirect detection 
is to use neutrino telescopes such as IceCube [26], ANTARES [27] and SuperKamiokande [28] 
to search for high-energy neutrinos produced by annihilation of WIMPs in the core of the 
Sun or Earth [29-33]. Due to their weak interactions with nuclei, such WIMPs would have 
scattered on solar nuclei and become gravitationally bound to the solar system, eventually 
returning to scatter repeatedly and settle (in the case of the Sun) to the solar core [34] (see 
also [35] for a review of this process and impacts of WIMPs on stellar structure, and [36-40] 
for additional corrections due to the influence of planets). 
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WIMPs arise in many proposed extensions to the SM. Standard neutrino telescope 
analyses [24-28, 41, 42] take an effective view of WIMP interactions, placing limits on WIMP- 
nucleon scattering cross-sections as a function of the WIMP mass, by assuming a certain 
annihilation cross-section and final state. These limits are difficult to translate into actual 
particle models, where the annihilation cross-section and branching fractions may take on a 
range of different values for any given WIMP mass and nuclear-scattering cross-section. To 
properly interpret limits on neutrino fluxes in terms of the parameters (defined at e.g. the 
Lagrangian level) of a theory for new physics, it becomes necessary to compare the observed 
neutrino flux with the predicted neutrino signal for each individual point in the parameter 
space of the theory. 

Having translated the flux limits into direct constraints on the parameters of a theory, 
it is then also possible to compare and combine the sensitivities of multiple experiments, 
even if they probe entirely different sectors of the theory (e.g. neutrino and accelerator 
searches). Within the realm of theories for new weak-scale physics, this 'global fit' approach 
has so far been applied mostly to supersymmetry (SUSY) [43-51], in the form of the minimal 
supersymmetric standard model (MSSM). Recent analyses have included new data from the 
Large Hadron Collider (LHC) [52-54], direct [55-58] and indirect detection [59-61]. Whilst 
great care needs to be taken over the details of the statistical methods employed [62-66], 
these analyses have proven a clear success, pointing the way to a future of closer comparison 
between astronomical and terrestrial experiments. 

To date global fits have not included neutrino telescope data. The ability of future 
incarnations of IceCube to detect WIMP annihilation in the constrained MSSM (CMSSM) 
has been studied [67-69] , based on the number of observed events only and simple estimates 
of the instrumental sensitivity. Other authors have looked at predicted neutrino rates in 2D 
slices through MSSM parameter spaces (e.g. [70-74]), or from sets of models in non-statistical 
random scans of various incarnations of the MSSM (e.g. [75-77]). 

Here we present a fast likelihood method for including full event-level data from neu- 
trino telescopes in global fits and related analyses. In particular, our formulation allows 
the directions and energy estimators associated with each event to be included in the final 
unbinned likelihood calculation, which can then be employed as a likelihood component in 
a global fit. The inclusion of spectral information in neutrino searches for WIMPs has been 
of particular interest recently [78, 79]. As a by-product of our approach, we also present 
a rigorous but simple exclusion measure for individual models, based just on the observed 
number of events in a neutrino telescope. We give full details of the input data required, and 
explicit examples using real data from the IceCube neutrino telescope. 

The IceCube data and simulations we employ are described in Section 2. These data 
have been made publicly available on the IceCube webserver [80] and in DarkSUSY 5.O.6. 1 The 
likelihood formalism, outlined in Section 3, is implemented and also available in DarkSUSY 
5.0.6. Our example global fits and MSSM scans are detailed in Section 4. We summarise our 
method and examples in Section 5. 

2 The IceCube Neutrino Telescope 
2.1 Description 

Completed on December 18 2010, the IceCube neutrino observatory [81] has 5160 digital 
optical modules (DOMs) installed on 86 strings between 1450 m and 2450 m below the surface 

www . darksusy. org 
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in the glacial ice at the South Pole. The detector consists of a hexagonal grid with horizontal 
spacing between strings of 125 m and a vertical spacing between DOMs of 17 m, leading to 
a total instrumented volume of 1 km 3 . Eight of the 86 strings were deployed in a more 
densely instrumented core in the middle of the array, with average inter-string separation of 
42 m and vertical DOM separation of 7 m. Together with the 12 adjacent standard IceCube 
strings, this forms the DeepCore subarray. DeepCore increases the sensitivity of IceCube 
at low energies, and substantially lowers the energy threshold. IceCube records Cerenkov 
light in the ice from relativistic charged particles created in neutrino interactions in or near 
the detector. By recording the arrival times and locations of these photons with DOMs, the 
direction and energy of the muon, and identity of its parent neutrino, can be reconstructed. 

2.2 Data samples 

For this paper, we use the analysis details of a search for WIMP dark matter annihilation in 
the Sun with the IceCube 22-string configuration [24, 80]. This data set has 104.3 days of live 
time, and was recorded between June 1 and September 23 2007. For the results in Sections 4.1 
and 4.2 we use the event-level data at final analysis level and signal simulations from [24]. As 
described in Ref. [24], final analysis level is reached through a series of increasingly stringent 
event selections that are applied to remove the cosmic-ray induced backgrounds. Events are 
selected which appear to originate from the direction of the Sun and which exhibit low energy 
signatures. 

The results discussed in Section 4.3 are based on a detailed study to determine the 
sensitivity of the 86-string detector to signals originating from dark matter annihilation in 
the centre of the Sun [41]. This study was performed as a full analysis in all details and gives 
a realistic expectation of the capabilities of IceCube to observe dark matter-induced sig- 
nals, given the state of data extraction, reconstruction, and signal discrimination techniques 
available at the time of the study. 

2.3 Signal and background simulation 

Previous work [24, 41] simulated solar WIMP signals using WimpSim [82, 83], which describes 
the annihilation of WIMPs inside the Sun. WimpSim simulates the production, interaction, 
oscillation and propagation of neutrinos from the core of the Sun to the detector, in a fully 
consistent three-flavour way. For the previous IceCube solar WIMP analyses that we consider 
(IceCube 22-string data analysis [24] and 86-string sensitivity analysis [41]), two annihilation 
channels were simulated: XiXi ~ * W + W~ (the 'hard' channel) and XiXi ~ * bb (the 'soft' 
channel). These two channels were chosen in an attempt to approximately cover the range 
of possible SUSY models, by assuming 100% branching into two channels with very different 
characteristics. For the 22-string analysis, neutralino masses m^o from 250 GeV to 5000 GeV 
were simulated. For the 86-string analysis, this range was extended down to m^o = 50 GeV. 
A similar exercise was performed for Kaluza-Klein WIMPs in models of Universal Extra 
Dimensions [25]. 

In contrast, the analysis method we describe here (Section 3) is developed to employ the 
specific annihilation spectrum appropriate for each model in a given theory of new physics. 
We use various parameterisations of the MSSM as examples, but the analysis technique we 
present is applicable to any theory containing a WIMP. To calculate neutrino annihilation 
spectra for individual SUSY models, we create yield tables from a series of WimpSim simula- 
tions, with neutralino masses from 3 GeV to 10 TeV and many different annihilation channels. 
For each SUSY model, we interpolate in these tables and add the contributions to the total 
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neutrino spectrum from each annihilation channel, according to the partial annihilation cross 
section for the given channel. We include more complicated annihilation channels involving 
model-dependent decays, such as those of Higgs bosons, by summing up the contributions 
from their decay in flight. This procedure is included in DarkSUSY [84], which we use to 
calculate the detailed neutrino spectra model by model. 

Background contributions for this analysis are muon events from single and coincident 
air showers, as well as atmospheric neutrinos. No dedicated background simulations are 
needed for the work we present here, as we estimate the expected background from scrambling 
real data at the final analysis level (detailed within Section 2.7). 

2.4 Effective area calculation 

We derived detector responses for the 22-string and 86-string IceCube configurations men- 
tioned in Section 2.1 from the sets of signal simulations used in each analysis [24, 41]. In 
each case, we calculated an effective area for detection of muon neutrinos by IceCube from 
the direction of the Sun as a function of neutrino energy. These effective areas correspond to 
averages over the austral winter. The 22-string area is available online [80]. For convenience, 
both datafiles are also redistributed in DarkSUSY 5.0.6. The effective area for muon neutrinos 
and muon anti-neutrinos are given separately. 

Detector systematics determined from the analyses are energy-dependent; the data files 
therefore include systematic uncertainties on the detector response within each energy bin, 
to which we assign a la confidence level. These uncertainties were determined within simu- 
lation studies, where identified sources of uncertainty, e.g. absolute DOM efficiency, photon 
propagation in ice or calibration constants, were individually varied within reasonable ranges 
of their original values. Similarly, the uncertainties arising from limited simulation statistics 
are also given for each energy bin of the effective areas, at the la confidence level. 

2.5 Angular response 

The point spread function (PSF) describes the uncertainty in the reconstructed arrival di- 
rections of incoming neutrinos. The full PSF is a 2D joint probability distribution for the 
angular separation between the true arrival direction and the reconstructed one, in each of 
two linearly-independent directions on the sky. For the small angles under consideration and 
suitably chosen axes, this is well approximated by a 2D Gaussian distribution of the form 



P(P u 2 ) = — !— exp 

Z7rCJlO"2 



el_9l 

2a\ 2a\ 



(2.1) 



where Q\ and #2 are the angular separation along each axis. Given that the Sun is a circularly 
symmetric source, we work exclusively with a reduction of the PSF into the 1-dimensional 
probability distribution for the absolute angular separation on the sky, <f> = \j6\-\-6\-> re- 
moving the azimuthal dependence. For an azimuthally-symmetric PSF, a = a\ = 02 and the 
ID PSF is 



p (4>) = -5 ex P 

a z 



2a 2 



(2.2) 



Over a given energy range, it is possible to accurately construct the ID PSF directly for 
neutrinos arriving from all directions, using IceCube data analysis signal simulations [24, 41] . 
In that case, one need not assume that the 2D PSF is even Gaussian, let alone azimuthally 
symmetric. However, the radial probability density function (PDF) of an azimuthally- 
symmetric 2D Gaussian is a reasonable approximation to the radial PDF of a moderately 
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non-symmetric 2D Gaussian, and has the added advantage that it can be parameterised in 
terms of a single a. As the 2D PSF in IceCube is not strongly asymmetric, we therefore 
model the ID PSF with Eq. 2.2, but extract the parameter a, which we refer to as the "mean 
angular error" , directly from the ID PSF constructed from IceCube signal simulations. This 
can be done by integrating to the containment angle implied by e.g. the mean ((</>) = yf-cr), 
second moment {{4> 2 ) = 2cr 2 ), or median (<7\/21n2; this is what we use here) of the (j) distri- 
bution in the signal simulations. As in previous work [24, 41], we consider a over the same 
energy bins as used for calculating the effective area. 

We associate angular uncertainties with real data events on an event-by-event basis, 
using the paraboloid method [85]. A paraboloid function is fitted to the muon track recon- 
struction likelihood function in the neighbourhood of the best fit. The resulting confidence 
ellipse on the sky is represented by the axes a\ and <72, which correspond to the standard de- 
viations of the likelihood function in each of two linearly-independent directions. The overall 
reconstructed likelihood track uncertainty, <r para (the "paraboloid sigma"), is calculated as 
the mean in quadrature of the uncertainties in the two axes, <7p ara = (o\ + <t|)/2. Good track 
fits generally result in a paraboloid that is narrow along both axes, and therefore have small 

<7para Values. 

The simulated distribution of reconstructed angles, used for determining the mean angu- 
lar error, contains contributions from both the intrinsic angular deviation of neutrino- induced 
muons from the incoming neutrino direction, and the error in the reconstructed muon di- 
rection itself. For actual data events, the former cannot be known except in an average 
sense, whereas the latter can be estimated for each event with the paraboloid sigma. For the 
22-string configuration, the reconstruction uncertainty dominates over the intrinsic angular 
deviation of neutrino- induced muons, so we are able to employ the paraboloid sigma £J para as 
an accurate approximation to the total angular error in Eq. 2.2, for individual events. For 
future IceCube configurations or other experiments where the muon angular reconstruction 
is of higher precision than the intrinsic muon angular spread, an additional component must 
be added to cx para for each event, in order to properly account for the difference between the 
incoming neutrino and neutrino-induced muon directions. 

2.6 Energy estimator 

In the broader context of IceCube data analysis, the study we present here deals with compar- 
atively low-energy neutrinos. The corresponding muon events are associated with minimum- 
ionising tracks, where stochastic losses, which increase with muon energy, are not dominant. 
Energy estimators, as generally used within IceCube analyses targeting higher muon ener- 
gies, exploit the energy-dependence of such losses in order to estimate the energy of the muon 
and the parent neutrino. The volume instrumented with the 22-string configuration is too 
small to use the reconstructed track length of either fully or partially contained muon tracks 
to determine the energy. Instead, we use the number of lit DOMs -/V c h an , a measure of the 
amount of recorded light per event, as a suitable energy estimator. We make no attempt 
to assign a specific energy and associated uncertainty to each event according to its iV c h an 
value. Instead, we calculate the expected distribution of observed -/V c h an values for a series of 
intervals in neutrino energy, and use these together with the predicted energy spectrum of 
the signal for a given SUSY model, to calculate the predicted distribution of -/V c h a n- Using 
an unbinned likelihood (Section 3), we then simply compare the observed value of -/V c h an 
for each event to this predicted distribution. We derived the probability distributions per 
z/-energy interval from high-statistics v simulations generated in the initial 22-string analysis 
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Figure 1. Predicted probability distributions of iV c han derived from high-statistics v simulations 
generated in the initial 22-string analysis [24]. Each distribution is constructed for neutrinos having 
energies in a specific logarithmic energy interval of width 0.2. The fitted Gaussian functions are to 
guide the eye only, and are not used in our calculations. 



[24]. These distributions are shown for the energy range relevant to our analysis in Fig. 1. 
Note that we do not use the fitted Gaussian functions in Fig. 1 for any purpose other than 
to guide the eye; we employ the actual distributions directly for our signal predictions and 
likelihood calculations. 

2.7 Background estimation 

We estimate the necessary background distributions from actual data. For the likelihood 
functions in Section 3, we need two background distributions: the angular distribution of 
background events d-Pec (</>') /d^>' — given as a function of 4>' , the angle between the recon- 
structed track direction and the Sun — and the distribution of iV c han due to background 
events, dP B G(AWi)/diVchan- 

The angular distribution of the dominant cosmic-ray shower background contribution 
is expected to be azimuthally independent, depending only upon the angle from the horizon 
(zenith angle). Thus, events incident in a horizontal band across the sky covering the same 
vertical range as the analysis (centred upon the solar zenith), but including all 360° around 
the horizon, are representative of the backgrounds in the small angular region around the 
Sun actually used in the analysis. This allows for an estimate of the necessary background 
distributions driven purely by observation. 2 The background is roughly constant per unit 
angular area in this band, with any mild zenith angle dependence softened by an averaging 

technically, the 360° band includes the direction of the Sun, so the events in this band may include some 
signal as well as background. However, the background overwhelmingly dominates any potential signal, so 
background estimates would be only negligibly affected. If this was a concern, a section of the band including 
the Sun can be removed and only an e.g. 300° band used. 
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over the zenith angle of the Sun during the austral winter. Hence ~j^, G ~ sin<^', though the 
actual d ^?, G we use is tabulated rather than fit to any functional form. 

2.8 Data format and availability 

The likelihood formalism described in the following Section, as implemented within Dark- 
SUSY 5.0.6, requires four simple text files containing IceCube data for any given analysis. 
The first data file gives the binned v and i> effective streets £is sl function of energy, along 
with the associated la statistical and systematic errors in each bin. This file also provides 
the mean angular error. The second data file contains all events that a) passed the event 
selection corresponding to the effective areas, and b) have reconstructed directions that fall 
within some specified cut cone around the solar position. The 22-string version of this file 
consists of 180 events within 10 degrees of the Sun. This file includes the observed value of 
A^ c han for each event, along with the reconstructed arrival angle relative to the Sun, and the 
corresponding paraboloid sigma. This file also includes the total live (i.e. exposure) time of 
observations towards the Sun. The third file gives the background distributions, with the 
angular distribution provided as a function of <fi' and the spectral distribution as a function 
of iVchan- The final file contains the iV c han response histograms presented in Fig. 1. 

We provide such files in DarkSUSY 5.0.6 for both the 22- and 86-string examples that 
we show in this paper. The 22-string files (containing real data), are constructed from 
the original IceCube analysis [24] and its associated data release [80]. The 86-string files 
(containing simulated data), are constructed from the IceCube sensitivity analysis [41] and 
the expected number of all-sky background events [86]. 3 It is anticipated that 79-string data 
will later also be made available on the web and in DarkSUSY, in connection with an analysis 
of those data using the techniques presented in this paper. 

3 Likelihood functions 

Using the data, responses and simulations described in Section 2, our goal is to evaluate 
to what degree IceCube observations support or constrain a theory of new physics ^, given 
some particular values ip of the theory's m free parameters. We refer to a specific new physics 
scenario $ as a 'theory' and a specific choice of its parameters ip as a 'model'. We may be 
interested in 

a) the absolute goodness of fit of the model with parameter set ip, or 

b) the goodness of fit of ip relative to other points in the parameter space of the theory. 

Case a) is model exclusion: we are interested in the maximal confidence level at which we 
can exclude the theory with specific parameter values ip as a true hypothesis. Case b) is 
most often parameter estimation: we are interested in determining how well ip fits the data, 
compared to the best-fit point ip within the parameter space of the theory. Assuming the 
overall theory to be correct, we use this information to determine a region of parameter space 
that includes the values of the parameters, with some specific level of confidence. 

The absolute goodness of fit for a model is described by its p- value, where the confidence 
level with which the model can be excluded is 1 — p. The p- value can be obtained from the 

3 In both cases the effective areas have been updated slightly with respect to the published values, to prop- 
erly account for neutrino-antineutrino asymmetries; these changes do not affect previous published IceCube 
limits. 
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likelihood function C(D\ip), where £(D\ip) dD describes the probability of observing data in 
the range [D, D + dD] if ^ is true and the true values of its parameters are ip. The relative 
goodness of fit can be obtained from the ratio 

AW) = (3-D 

Here £ m ax is the maximum likelihood, found at the best-fit point ip. Calculation of p or A re- 
quires knowledge of the distribution of C(D\tp) for repeated experiments (i.e. the distribution 
in data space). The simplest approach is to assume that —2 In A follows a \ 2 distribution 
with m degrees of freedom (as predicted in the limit of infinite data by Wilks' Theorem [87]), 
whilst the most rigorous approach is to explicitly construct the distribution by brute force 
simulation (the Feldman-Cousins approach [88]). 

For both model exclusion and parameter estimation, we require an expression for 
£(D\ip). Here we develop a likelihood function for IceCube data, based on the total number 
of observed events and the individual properties of each event. We go on to develop a p- value 
for model exclusion based on our likelihood construction. Although we work exclusively at 
neutrino level in this paper (cf. Section 2.2), the treatment in this Section should be equally 
valid at muon level. 

3.1 General unbinned likelihood 

Consider a set of n to t events observed in IceCube. Denote the true energy of the ith event as 
E{, and its arrival angle relative to the position of the Sun on the sky as fa. The ith event 
will be reconstructed with arrival angle <\>\, and will cause a number of DOMs iVj to fire. The 
unbinned likelihood for these data is 

'^tot pTT poo dP 

£unbin = £(ntot|0tot) TT / / Q{N i ,(f>' l \E i ,cP i )———(E h &,7P)dE i d<f )i . (3.2) 

fJ^Jo Jo dEidcpi 

Here Q(Ni,(j/ i \Ei,(j)i) is the probability per unit angle and iVchan (i.e. probability density) for 
observing N{ and <^ for the ith event when the true values of the energy and arrival angle 
are Ei and <j>{. The a priori probability density for the ith event to arrive with energy E± 
from angle fa is given by , . (Ei, V0> this is a prediction of the model i/j- 

3.2 Number likelihood 

Neglecting systematic issues, which we address below, the prefactor £(ntot|#tot) in Eq. 3.2 is 
the standard Poissonian likelihood, equal to the probability of observing ntot events produced 
by a Poisson process with mean 9 tot 

antat p — #tot 

£(ntot|0tot) = tot , • (3.3) 

"-tot! 

The expectation value for the number of events is #tot = ^tot(V ; )) the total number of events 
predicted by model ip. This is given by 



'tot 



(V)=0BG + 0sty), (3-4) 



the sum of the predicted number of background events #bg and signal events #s(V0- Whilst 
the predicted number of signal events depends on the model, the background does not. 
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Eq. 3.3 accounts for statistical fluctuations in the number of observed events for a 
precisely- known mean 6 to t- However, the predicted mean itself may contain a systematic 
error, due to e.g. an error in the estimate of the effective area of the instrument. Taking e to be 
the (unknown) ratio of the true expected signal contribution 0$ true to the nominally- predicted 
signal contribution 9s, we have #s,true = e#s- This means that the relative (fractional) error 
on #s is e — 1- To account for this potential systematic error in the number of predicted events, 
we marginalise over a probability distribution for the relative error e — 1 in a semi-Bayesian 
manner [89] . Assuming a log-normal probability distribution for e — 1 , the likelihood becomes 



<£num ("tot I #tot 1 



1 



2-K(J f 



"-tot! 



exp 



1 /lne 

2 1^7 



de. (3.5) 



Here a e is the fractional systematic error in the predicted number of events. We take this term 
to be the sum in quadrature of the fractional uncertainty of the IceCube effective area, and 
a theoretical error r. The fractional error in the effective area is itself the sum in quadrature 
of the fractional statistical and systematic errors determined in deriving the effective area 
(both have the character of a systematic for the purposes of the likelihood calculation, if not 
for the effective area calculation itself). For r we use a conservative error of 5%, designed to 
account for neglected higher-order corrections and possible accumulated numerical round-off 
errors. 

We apply the systematic uncertainty only to the signal prediction in Eq. 3.5, as neither 
an uncertainty in the signal prediction from theory, nor an error in the estimation of the 
effective area, would impact the number of predicted background events. This is because 
#BG is based on the total number of observed events away from the Sun. Eq. 3.5 gives the 
complete contribution of the number of events to the total likelihood; we hence refer to this 
as our adopted 'number likelihood'. 

We could have adopted a Gaussian probability distribution for e — 1 in Eq. 3.5, but this 
leads to a non-zero probability that e = 0, which is unrealistic. We have therefore adopted 
a log- normal distribution for e in Eq. 3.5 and throughout this paper, although for small o~ t 
the difference with the Gaussian case is minimal. As the Gaussian calculation is significantly 
faster in some cases, we provide both Gaussian and log-normal routines in DarkSUSY 5.0.6. 

Eq. 3.5 assumes a single, constant systematic error on the effective area. For real 
detectors, this error is typically highly energy-dependent. In general there is no consistent way 
to allow the prior distribution for e in Eq. 3.5 to vary with energy (or any other observable) 
in an unbinned likelihood calculation. One option is to model the variation of the systematic 
error with energy using some parametric form, and then marginalise over the parameters; 
this however requires assumptions as to the functional form and the degree to which the 
size of the systematic is correlated at different energies. Another solution is to bin events 
into broad 'energy' ranges according to their observed Af c h an values, using a Bayesian sorting 
technique that assumes a prior for the source spectral shape. Each bin is then assigned a 
different systematic error for the effective area, based on the range of neutrino energies it 
covers, and an unbinned analysis is performed on the events within each bin. Ideally, one 
keeps the number of bins in such a setup to a minimum, to minimise prior-dependence in 
the final results. A third, more conservative option, is simply to take a e in Eq. 3.5 to be 
the largest percentage error seen on the effective area at any energy. We follow the third 
strategy; we also attempted an analysis using the second, but found that the additional noise 
introduced by the need to sort events into the different bins grossly outweighed any advantage 
gleaned from having a more accurate effective area at higher neutrino energies. 
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3.3 Spectral and angular likelihoods 

In general Q{Ni,fa i \Ei,(pi) in Eq. 3.2 is a function of the spectral and angular resolution of 
the detector 

Q{N h cj ) ' i \E i ,fa^) = E^N^PSF^fa), (3.6) 

where -Edisp is the distribution of A^an for incoming neutrinos of energy Ei (the energy 
dispersion of the instrument), and PSF is the detector point-spread function. Here we have 
assumed that E'disp has no angular dependence, and that the energy-dependence of the PSF 
can be neglected at this stage. These assumptions are certainly not true in general, given 
the large width of E^ sp and the strong energy-dependence of the PSF. Here we apply a hard 
angular cut around the solar position, reducing the possible range of fa, and therefore the 
range of variation of -E"disp with angle. Although this does not make E^ sp entirely isotropic 
even within the analysis cone, our use of such a cone makes neglecting this angular variation 
reasonable. The second approximation is well justified because our per-event analysis already 
takes the energy-dependence of the PSF into account implicitly, because we apply a unique 
PSF for each event using the paraboloid sigma (cf. Section 2.5). 

Similarly, d& d</>. {F%, fa,/fa) in Eq. 3.2 depends on the true energy spectrum ^ and 
spatial distribution ^ of the signal and background 



dP dP dP 

-{E l ,fa) = —{E l ,i ) ) — {fa^), (3.7) 



d£i dfa dEi dfa 



with 



g^^) = /BC^f(^ + /sg(E^) (3-8) 

^{fa^) = fBG d ^{fa) + fs^{fa,^- (3.9) 
dfa dfa dfa 

Here d ^ G , , etc refer to the signal and background contributions to the overall predicted 
spectral and angular distributions, and 

A = AW = F^. /bg = /bgW = /jtt (3.10) 

are the signal and background fractions, with 6s and 6bg the total predicted number of signal 
and background events within the analysis cone. Like #tot and 6s, the fractions fy and /bg 
depend upon the model ip by definition, but we typically do not write this dependence out 
explicitly. 4 

In Eqs. 3.7-3.9 we have again dropped any explicit energy dependence in the angular 
part, or angular dependence in the spectral part. This is perfectly well justified for the 
background component, as the arrival directions and A^han values (and therefore by inference, 
energies) of background events are observed to be essentially uncorrelated. For the signal 
component alone, the assumption that the energy spectrum does not vary with arrival angle 
is well justified because the IceCube PSF is far larger than the the angular extent of the 



4 Note that although not all the actual background is genuinely due to neutrinos (there is a large atmospheric 
muon component), we can simply interpret the parts of Eqs. 3.7-3.9 due to non-neutrino events as effective 
neutrino angles and energies corresponding to those tracks. As we only ever deal with the observed background 
distributions, not the predicted backgrounds (cf. Eq. 3.f7), this is entirely valid. 
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region where high-energy neutrinos are produced in the Sun, regardless of the energies of 
the neutrinos. That the angular extent of the signal does not depend on the energy of the 
neutrinos from the solar core, nor therefore the model ip, follows on exactly the same grounds. 
(We can therefore model the angular distribution of the signal as simply a delta function at 
the solar position.) 

However, these assumptions do not rigorously hold for the combined signal and back- 
ground prediction, as the spatial and spectral characteristics of the signal and background 
differ; neutrinos coming from the Sun should presumably exhibit a more signal-like spectrum 
on average than those arriving from elsewhere on the sky. We deal with this complication 
implicitly by way of our finite analysis cone. We evaluate the signal spectral distribution 
in Eq. 3.8 at precisely the solar position, and then weight its contribution against that of 
the background by way of the total predicted signal and background fractions fs and /bg- 
These fractions are the integrated predictions over the full analysis cone, meaning that in 
Eq. 3.8 we essentially derive a mean predicted energy spectrum over the full cone. Similarly, 
we implicitly account for the energy-dependence of the angular distribution by weighting the 
signal and background angular distributions by the total signal and background fractions. 

Collecting all the terms from Eq. 3.2 where we have retained an explicit dependence on 
the energy of an incoming neutrino or the number of DOMs it triggers in the detector, we see 
that the contribution to the overall likelihood from the spectral properties (i.e. iVj) of event 
i is 

poo d p 

C spec ,i( N i\^) = / E disp (N i \E i ) — (E i ,2l;)dE t (3.11) 
Jo d-frj 

/BG^disp(A r ^i)^^(^) + /s^dis P (iV i |^)^(^ 5 V') dE t . (3.12) 

The true energy spectrum of the atmospheric background d ^ G would need to be estimated 
empirically using the observed background flux away from the Sun. The distribution actually 
observed by IceCube (Section 2.7) is the true background distribution d ^ G convolved with 
-^disp! this is simply the first term in Eq. 3.12, so 

dPnr f°° dPs 

£speci(W) = fBG-rj^(Ni) + fs / £dis P (^)-7#(£i, VO d^. (3.13) 

Collecting the remaining terms in Eq. 3.2 with an explicit direction-dependence, the 
angular likelihood for event i is 

f n dP 

^MW)= I PSFttM — ttu^dfr (3.14) 

fBG^^^PSF^) + fs^^^PSFi^dcP,. (3.15) 
d(pi d(pi 

As for the spectral distribution of the background, we estimate the angular dependence of 
the background empirically, based on observations away from the Sun (cf. Section 2.7). This 
provides a direct estimate of the first term in Eq. 3.15. 

The predicted angular distribution of signal events depends on the actual angular 
dependence of the signal. Given that the Sun is a point source to IceCube, we can simply write 



OO 



- 11 - 



dPs/d(j)i = S((j>i). We can therefore easily solve the integral in the second term in Eq. 3.15. 
Combining this with the known angular distribution of background events, Eq. 3.15 becomes 

£.n S M\i>) = fsPSF^i = 0) + / BG ^(^). (3.17) 

PSF(4>' i \(j) i ) in Eqs. 3.6, 3.14, 3.15 and 3.17 refers to the unique ID reduced PSF for 
each event i, as determined by its paraboloid sigma Oi. With a full PSF Gaussian in 2D (cf. 
Section 2.5) and an analysis cone cut angle 0cut) the reduced PSF is given by 



PSF{(j)' i \(l) i ) = C(0i, ai, <#.ut) ' 2 ex P 



1 



2 V ai 



where 



C(4>i,ai,(j)' CVLt ) =l-exp 



1 



(-'cut 



2 \ Gi 



(3.18) 



(3.19) 



is an angular correction factor that ensures the PSF is unitary within the analysis cone (i.e. 
the integrated probability for any given event to have come from anywhere at all is 1). 

3.4 Total likelihood function 

To arrive at the final composite unbinned likelihood, we simply substitute Eqs. 3.5, 3.13 and 
3.17 into Eq. 3.2, giving 

"tot 

^total^tot; ^IVO = -£num(^tot 

i=l 

where S = {iVj, ^}i=i..ntot * s the set of all reconstructed event arrival directions and -/V c h an 
values. This is the likelihood function we employ for all the examples of parameter estimation 
in this paper. 

3.5 Predicted distributions and event counts 

So far we have focused on the construction of the likelihood function, assuming that one 
knows how to calculate the number of events and their distribution in energy predicted by a 
model The signal contains contributions from both neutrinos and anti-neutrinos, 

e s y>) = 6sA4) + 0sA4)- ( 3 - 21 ) 

The energy distribution of the predicted signal is the normalised sum of the neutrino and 
anti- neutrino flux spectra indicated by model ip, each weighted by the respective effective 

For an extended source, one also needs to take into account any anisotropy in the exposure time t exp (<f>) of 
the detector, defined as the total time over which the detector is sensitive to neutrinos arriving from true angle 
4>. The effective exposure may vary across the sky due to the seasonal variation of the detector orientation, 
and choice of analysis cuts. In this case, the predicted distribution would be 

xrOW = ^MWfc) / r ^(0)i cxp (0)d0 , (3.i6) 

where dPs,i/d<^i is the intrinsic source distribution, and the denominator simply ensures correct normalisation. 
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area of the detector (A u or Ap). The sum is further weighted by an overall angular loss factor 
L, giving 



~dE ' = 0s L ^'^ c ' 



(3.22) 



Here d ^ l/ VO an d (-£?, ip) are the respective differential neutrino and anti-neutrino 
fluxes incident on the detector (due to neutrino production in the Sun), with units of events 
per unit area, energy and time. The prefactor t exp ((p = 0)/9$ tl , ensures that Eq. 3.22 is 
indeed normalised. Both A and L are explicitly energy-dependent; A has units of area and 
L is dimensionless. 

The angular loss factor L(E, (fi' cut ) accounts for events that originate from the Sun, but 
are so badly reconstructed that they fall outside the analysis cut cone ((// > (/>' cut ). We use 
the (ID, reduced) mean PSF of the instrument to determine L, as defined by the energy- 
dependent mean angular error &{E) (cf. Section 2.5). This gives 



L(E, 0' cut ) = 1 - exp 



1 ' ' 'cm 



/ \ 2 



2 \a{E) 



(3.23) 



We can simply read off the definition of 6s iU from Eq. 3.22 by virtue of it being a 
normalised distribution 



= 0) J™ L(E)A V (E)^^(E, V0d£. (3.24) 



The corresponding expression for 6$ t p is Eq. 3.24 with v — > v. This integral has no cutoff, as 
lim^o A(E) = lim^oo A(E) = 0. ' 

The total number of predicted background events in the analysis cone (#bg) can be 
determined by simply rescaling the background rate observed away from the Sun, to the 
desired exposure time and analysis cone sky fraction. 

3.6 p-value for model exclusion 

Unfortunately, there is no established method for determining goodness-of-fit using unbinned 
maximum likelihood estimators, short of a full Neyman construction using Monte Carlo 
simulations based on Eq. 3.20. Whilst this is possible for simple toy models, it is not a 
computationally feasible option for global fits in virtually any UV-complete theory; lengthy 
renormalisation-group and relic density calculations make it unrealistic for analyses of super- 
symmetry, for example. 

A more realistic option is to do model exclusion based entirely on the number of observed 
events, using Eq. 3.5. This has the obvious disadvantage of discarding the spectral and 
angular information that go into Eq. 3.20, reducing our overall ability to exclude models 
compared to a full Neyman construction. However, this strategy has the distinct advantage 
that we know exactly how Eq. 3.5 is distributed in data space, as it simply describes a 
Poisson process. This knowledge is what allows us to calculate an absolute p- value for any 
model, completely independent of the value of the likelihood function elsewhere in a theory's 
parameter space, and exclude that model with confidence 1 — p. 

We begin by identifying each ifj as a specific signal hypothesis, and identifying ifj plus 
the background distribution inferred from observations away from the Sun as a specific sig- 
nal+background hypothesis s + b. The p-value for a hypothesis is defined [90, 91] as the 
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probability, in a repeat of the chosen experiment, of obtaining a test statistic at least as 
extreme as the one actually observed, if the hypothesis is true. The test statistic with which 
we are probing the s + b hypothesis is the likelihood ratio, 

x _ £num(rctot|fls + Abg) ,g 25) 

Amm(rttot|0BG) 

X is the ratio of the number likelihood (Eq. 3.5) for the signal+background prediction to the 
number likelihood for a background-only prediction. A more extreme result is one that gives 
a lower likelihood ratio (i.e. is less probable if s + b is true). We do not probe the hypothesis 
ip directly with Eqs. 3.5 and 3.25, as ip alone does not give a specific prediction for the total 
number of observed events. 

Referring to Eqs. 3.3, 3.5 and 3.25 we see that X is in fact a monotonically increasing 
function of the number of observed events ntot- This means that the p- value for a hypothesis, 
as tested using the number likelihood ratio X, is the total probability of observing ntot 
or less events if the hypothesis is true. This is precisely the sum of the likelihood over all 
possible numbers of observed events less than or equal to n to t- For the signal plus background 
hypothesis, this is 

"tot 

Ps+b = ^ ^num [n|0tot (VO] ' ( 3 - 26 ) 
n=0 

and for the background-only hypothesis 

"tot 

Pb = > / nuln (n|0 B G) • (3.27) 
n=0 

Note that there is a small conceptual complication in the case of the background-only hypoth- 
esis: direct substitution of 0s = into Eq. 3.25 causes the distribution of the test statistic 
to collapse to a delta function at X = 1. Eq. 3.27 is therefore strictly valid only when 
Eq. 3.25 is taken in the limit 0s ~~ >■ 0, rather than evaluated at exactly 0s = 0. For further 
explanation, we refer the interested reader to the original description [92, 93] of the CL S 
technique (of which our p- value is essentially an example), and the conditioned confidence 
interval technique from which the CL S method was generalised [94]. 

For a true Poisson likelihood given by Eq. 3.3, the sums in Eqs. 3.26 and 3.27 can be 
performed relatively painlessly, but for the "smeared" Poisson likelihood given by Eq. 3.5 (the 
form we actually use 6 ), the sum is over terms each requiring a numerical integration, a po- 
tentially slow process. We provide optimised routines for handling the numerical integration 
and summing that ensure calculations of the p- values are relatively quick. 

In a classical sense p s +b can be regarded as the p-value for testing the hypothesis ip 
when the background hypothesis is known to be correct. However, it is sensitive to statistical 
fluctuations in the actual observed background; a downward fluctuation in the background 
can cause ip to be excluded with an unreasonably high confidence level, for example. To 
correct for such problems, we use the modified p- value as our estimate for the hypothesis ip 
alone, defined as [90] 

Jty = ^. (3.28) 
Pb 



6 Strictly speaking, it is Eq. 3.3 that we know the distribution of exactly, not Eq. 3.5. Given that Eq. 3.5 
is simply a smeared form of Eq. 3.3 however, the two distributions should be very similar. 



- 14 - 



4 Example SUSY scans and global fits 

As concrete examples of our likelihood formalism in action, we perform three different anal- 
yses of the impact of IceCube data on theories of new physics. We examine two variants of 
the minimal supersymmetric standard model (MSSM): the constrained MSSM (CMSSM), 
and a 7-parameter, low-energy phenomenological parameterisation (the MSSM-7). In the 
first analysis, we consider the impact of including real IceCube 22-string data on a CMSSM 
global fit, using the full likelihood function (Eq. 3.20). In this first analysis, as a rough 
proxy for the full detector we also rescale the effective area of the 22-string analysis, and 
consider the resulting impact on the CMSSM. In the second analysis, we assess the ability 
of our method to recover a hypothetical WIMP annihilation signal from the Sun in terms 
of the CMSSM, using the same rescaled 22-string detector. In the third analysis, we give 
an example of a simple model exclusion exercise in the MSSM-7, using Eq. 3.28 with the 
86-string detector simulation. 

4.1 CMSSM global fit with 22-string data 

The CMSSM is a high-energy parameterisation of the MSSM. It is defined in terms of the 
sign of the Higgs mixing parameter fj,, the ratio of up-type to down-type Higgs VEVs tan/3, 
the universal trilinear coupling Aq and the universal mass parameters mo and rn\/2- tan/3 is 
defined at the electroweak scale, whereas Aq, mo and m^j are defined at the scale of grand 
unification. The WIMP mass, given by the mass of the lightest neutralino m v o, is a derived 

Al 

quantity. We refer the reader to Ref. [95] for full details. 

In order to include IceCube 22-string data in a CMSSM global fit, we coupled Dark- 
SUSY 5.0.6 to SuperBayeS 1.5.1 [46, 47, 67] and included Eq. 3.20 as an additional likelihood 
component in the global fit. The total likelihood function we employ for scans in this section 
includes all components available in the public release of SuperBayeS 1.5.1: LEP limits on 
sparticle and Higgs masses, the anomalous magnetic moment of the muon, limits on rare 
processes like b — > sj, additional S-physics observables, the dark matter relic density, and 
measurements of SM nuisance parameters (refer to [47, 67] for details). We do not include 
more recent constraints from the LHC [51-54, 96, 97], XENON-100 [57, 58, 61] or indirect 
detection [59-61], as they are not necessary for our simple illustrative purposes here. Future 
analyses using the likelihood formalism we describe, where the focus is on physical interpre- 
tation rather than methods, should include such data. 

We performed a scan of the CMSSM parameter space using MultiNest [98, 99] with 4000 
live points and a tolerance of 0.5. These are settings appropriate for mapping the posterior 
PDF in a CMSSM scan, but not the profile likelihood [62, 63]. The likelihood formalism we 
present here is of course equally well adapted for use in frequentist as Bayesian analyses, but 
obtaining correctly converged profile likelihoods takes roughly an order of magnitude more 
computing power than posterior PDFs. We thus rely almost exclusively on posterior maps 
in this paper, as we are essentially only interested in the results of the scans for the sake of 
illustration, so the choice of posterior PDF or profile likelihood is fairly arbitrary. One must 
also be concerned with correct coverage in presenting profile likelihood results [64-66], an 
additional unnecessary complication for our current purposes. 

An important point however is that the full IceCube likelihood calculation takes a 
fraction of a second in nearly all cases, so it is not a significant bottleneck in most scans 
when compared to e.g. a relic density calculation. 
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Figure 2. Joint 2D posterior probability distributions from a CMSSM global fit including IceCube 
22-string event data. Contours indicate 68.3% (la) and 95.4% (2a) credible regions. Shading and 
black contours give the result for the fit with IceCube data, whilst grey contours correspond to an 
identical global fit, but with no IceCube data included. Existing 22-string data has a minimal impact 
on the parameter space of the CMSSM. 

We parsed and plotted the results of our SuperBayeS runs using pippi [100]. The parsing 
step requires binning samples and marginalising posterior PDFs or profiling the likelihood 
function. In most cases, we sorted samples into 100 bins in each parameter /observable, 
and interpolated between those bins for display, at a resolution of 500 points per parame- 
ter/observable. Plots of the predicted number of signal events in IceCube are the exception: 
in this case we used 130 bins, and interpolated to a resolution of 3000 points per observable 
direction, to better resolve the region of interest. 

Our first global fit employed the observed events, detector response and effective area 
of the published 22-string (IC22) analysis [24]. For this analysis we used an angular cut of 
4»cut = 10° around the solar position. The resulting marginalised posterior PDFs are are 
shown in colour in Fig. 2, along with contours indicating 68.3% (la) and 95.4% (2a) credible 
intervals. Here we show results in the standard mo, m^ parameter plane, as well as in terms 
of the neutralino mass - nuclear scattering cross-section planes, for both spin-dependent and 
spin-independent interactions. We also plot corresponding credible contours in grey for an 
identical scan run without the inclusion of any IceCube data. 

As expected from the limits produced in the original 22-string analysis [24], this data 
has very little impact upon the preferred regions in a CMSSM fit; the stau co-annihilation 
region (the small region at low mo in the leftmost panel and the curved, lower cross-section 
regions in the other two panels) is unaffected, and any impact on the focus point region (the 
larger region at high mo in the leftmost panel and the nearly horizontal, high cross-section 
bands in the other panels) is difficult to see. On very close comparison of the black and 
grey contours in the right panel of Fig. 2 (somewhat easier if one actually refers to the grey 
contours in the right panel of the next figure, Fig. 3, which are identical to the grey contours 
plotted in Fig. 2), a very small impact can perhaps be seen. A tiny part of the focus point 
region, at large spin-dependent scattering cross-section and m v o ~ 200 GeV, appears to have 
been disfavoured by the IC22 likelihood at the level of moving it from the la credible region 
to the 2a region. Whether this is simply scanning noise is difficult to say, but the fact that 
the region in question has a very high posterior PDF in the scan without IceCube data, and 
is therefore very well-sampled, would seem to argue against such an interpretation. 

To illustrate the approximate effect on the CMSSM of an eventual non-detection of 
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Figure 3. As per Fig. 2, but assuming that the detector effective area is boosted by a factor of 100, 
in order to make a rough estimate of the impact of the complete IceCube detector on the CMSSM. 
Grey contours again correspond to the fit with no IceCube data. Upcoming searches for dark matter 
with IceCube will robustly exclude the majority of the focus point region of the CMSSM. 

neutrinos from WIMP annihilation in the Sun by the full IceCube detector, using our like- 
lihood formalism, we performed a second global fit with a rescaled IceCube effective area. 
We multiplied the IC22 effective area by a factor of 100, and kept all other aspects of the 
detector as in the 22-string analysis (angular errors, event sample, backgrounds and spectral 
response); we refer to this as the 'IC22 x 100' configuration. Although we employ the actual 
simulated 86-string analysis [41] later for model exclusion, using it for a study such as we 
describe in this section is not possible, as the 86-string analysis does not contain the requisite 
A'chan information to include spectral information in the likelihood function. 

The results of the IC22 x 100 global fit are shown in Fig. 3. Grey contours again 
refer to an identical scan performed without the inclusion of any IceCube data. As is clear 
from these results, the sensitivity of something resembling the full IceCube detector to both 
spin-independent and spin-dependent WIMP-nucleon interactions should place very strong 
constraints on the CMSSM, all but ruling out the majority of the focus-point region. Whilst 
this has been shown to already be the case when XENON- 100 data is included in the global 
fit [57, 58], those constraints are based only on spin-independent scattering, and are therefore 
particularly sensitive to the adopted prior on the hadronic matrix elements [101]. If the value 
of the pion nuclear sigma term is taken from lattice calculations instead of from experiment, 
XENON-100 provides almost no constraint on the focus point region. IceCube should provide 
a more complete exclusion of most of the focus-point region than XENON-100, because 
both spin-independent and spin-dependent couplings are expected to contribute to the solar 
capture rate, and spin-dependent scattering is far less sensitive to hadronic uncertainties than 
spin-independent scattering. The full IceCube result will therefore constitute an important 
independent verification of the XENON-100 exclusion; if IceCube sees a signal in its 86-string 
configuration, this will be a strong indication that the CMSSM is not responsible for dark 
matter, or that there is an error in the experimental determination of the pion nuclear sigma 
term. Even if there is no signal, IceCube will have a major impact on the parameter space of 
less constrained versions of supersymmetry, where the spin-independent and spin-dependent 
nuclear scattering cross-sections are not so tightly coupled as in the CMSSM. 

The appearance of the narrow region at low TO1/2 an d m x o (the well-known CMSSM 
light Higgs funnel region) only in the IC22 x 100 scan is not surprising, and can be understood 
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via a combination of prior and scanning effects. This region 'exists' as a locale of reasonable 
fits in all scans, but only appears in the IC22 x 100 posterior plots because in this case, 
given the linear prior we have employed on and mo, removing most of the focus point 
increases the relative contribution of the funnel to the final model evidence. These nuances 
are not directly related to the inclusion of IceCube data, and have been examined extensively 
already in the literature [47, 62, 63]. 

4.2 CMSSM recovery validation with mock IceCube data 

To explicitly test the validity and performance of the likelihood construction outlined in 
Sec. 3, we performed a signal recovery within the CMSSM, using the IC22 x 100 detector 
configuration. We simulated a WIMP annihilation signal in the Sun expected to produce, on 
average, 60 signal events within the 10 degree angular cut. We did this as a single realisation of 
the angular and spectral distribution expected from a 500 GeV WIMP annihilating exclusively 
into W + W~ pairs, taking into account the angular resolution of the 22-string detector and the 
expected distribution of iV c han values from such a signal. We added this signal to a realisation 
of the expected background, which we obtained by bootstrap Monte Carlo simulation using 
the actual all-sky events observed in the 22-string analysis (as described in Sec. 2.7). The 
particular realisation we employed contained 53 signal and 164 background events within 10 
degrees of the solar position. 

We then ran a full parameter scan based on the simulated event list, using DarkSUSY 
5.0.6 and our modified version of SuperBayeS. For this recovery we included only the sim- 
ulated IceCube data in the likelihood function (i.e. no relic density or other experimental 
constraints) , as we are specifically interested in the ability of IceCube to pin down parame- 
ters of a theory of new physics using our likelihood formalism. 

In Fig. 4 we show the results of the recovery in terms of various two-dimensional 
marginalised posterior PDFs. Here we give not only joint distributions for mo and mj/2, 
and for m x o and the two nuclear-scattering cross-sections, but also for m x o and the velocity- 
averaged annihilation cross-section {av)o, and m x o and the number of predicted signal events 
in IC22 x 100. We also give in Fig. 5 the ID marginalised posterior probability distributions 
of the WIMP mass m v o, and the three cross-sections. For comparison, we give the results of 
a similar recovery not including the spectral component of the likelihood in grey. 

Figs. 4 and 5 show that that our technique accurately recovers both the WIMP mass and 
expected number of events, and results in rather accurate measurements of the various cross- 
sections. Fig. 4a demonstrates that in this particular example, the fit hones in quite effectively 
on a part of the focus point region, corresponding to quite small regions in the mass-cross- 
section planes (Figs. 4b,c,d). The performance in the cross-section directions (Figs. 4b,c,d 
and 5b,c,d) is particularly impressive given that our benchmark point does not specify any 
cross-sections, only an expected number of events. A number of different combinations of 
cross-sections result in an appropriate number of events, as capture and annihilation need 
not be in equilibrium in the Sun in our calculation, and the relative contributions of spin- 
independent and spin-dependent scattering to the total capture rate is not specified. 

Although the benchmark point lies just outside the edge of the la posterior credible 
region in m x o and the predicted number of signal events (Fig. 4e), this is entirely consistent 
with realisation noise and the impact of priors. The profiled likelihood function (Fig. 4f) in 
this plane easily encompasses the true value within its la confidence interval, indicating that 
the slight preference towards lower masses in the posterior PDF is a prior-driven effect. As 
expected, Fig. 4f is somewhat noisy because the scan was not optimised for profile likelihoods, 
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Figure 4. 2D posterior probability distributions (a - e) and one comparison 2D profile likelihood (f), 
from a CMSSM signal recovery based on IceCube 22-string event data only. We have again assumed 
that the detector effective area is boosted by a factor of 100 relative to the regular 22-string analysis, 
to approximate the final detector configuration of IceCube. Here we have injected a simulated signal 
corresponding to a 500 GeV WIMP annihilating into W + W~ , at a rate that would give 60 signal events 
on average. Contours indicate 68.3% (ler) and 95.4% (2a) credible/confidence regions. Shading and 
black contours indicate the recovery using the full likelihood (Eq. 3.20), containing number, angular 
and spectral information. Grey contours give the recovery achieved with only number and angular 
information. The recovery is quite tight, especially with the inclusion of spectral information. 
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Figure 5. ID posterior probability distributions for the mass of the lightest neutralino and its most 
important cross-sections, from a CMSSM signal recovery. This is the same recovery shown in Fig. 4, 
based on IceCube 22-string event data only, but assuming a factor of 100 boost in the effective area in 
order to approximate the final detector configuration. Blue curves give the probability distributions 
obtained using the full likelihood (Eq. 3.20). Grey curves show the distributions achieved without the 
use of spectral information. The use of spectral information significantly improves the accuracy with 
which the mass and cross-sections can be determined. 

but the convergence towards a roughly Gaussian likelihood in this plane is obvious. Recall 
also that the benchmark point refers to the expected number of total signal events (60), but 
that the random realisation that we happen to employ contains 53 signal events inside the 
analysis cone. It is therefore not at all surprising to find that the best-fit and posterior mean 
values returned by the scan lie in the range 45-50, given the size of the credible/confidence 
regions. 

The grey contours in Fig. 4, and corresponding grey curves in Fig. 5, illustrate the dra- 
matic improvement in the mass recovery achieved by including even a crude energy estimator 
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Figure 6. A model exclusion analysis in the MSSM-7, based on Eq. 3.28 and the consideration of all 
events seen to arrive from a direction within 20 degrees of the solar position. Here we use a simulated 
86-string IceCube configuration, and simulated data consistent with background. A number of models 
will be directly excludable by the full IceCube detector, without any recourse to the rest of the model 
parameter space, nor any spectral or angular information beyond the number of events within the 
angular cut cone. The strongest exclusions will be for models with mixed gaugino/Higgsino dark 
matter, giving rise to large spin-dependent nuclear scattering rates. In particular, many of these 
models have spin-independent cross-sections too low to be detected by large terrestrial experiments 
like XENON-100. 

like iV c han- We see that the full likelihood consistently disfavours low masses compared to the 
recovery without spectral information. This results in smaller credible regions in all planes 
of Fig. 4, particularly those involving cross-sections, and tighter credible intervals in Fig. 5. 
Unsurprisingly, excluding the angular part of the likelihood also has a dramatic effect; using 
just the number likelihood results in a 2a credible region that covers most of the parameter 
space. 

4.3 MSSM-7 model exclusion with 86-string simulation 

We have also performed scans in the more phenomenological MSSM-7 with seven free param- 
eters at the electroweak scale: the Higgsino mass parameter /z, the gaugino mass parameter 
M2, the ratio of the Higgs vacuum expectation values tan/3, a common sfermion mass pa- 
rameter mo, the CP-odd Higgs boson mass txia and the trilinear couplings for the third 
generation squarks Af, and At- See e.g. [102] for a more detailed description of this model. 

The database we use here is made from several scans focusing on different regions of 
the parameter space. The core of the database consists of general random scans allowing the 
mass parameters to take on values up to lOTeV, At and Af, in the range between — 3mo and 
3mo and tan/3 between 1.2 and 60. We have also performed importance sampling scans with 
ADSCAN [103]. ADSCAN is setup to use an adaptive integration routine (VEGAS) based 
on importance sampling. We have then chosen to run ADSCAN with different "importance 
functions" G focusing on different phenomenologically interesting areas in the parameter 
space. The importance function we have mostly used is 

G {n x h 2 ) = exp 

with fiwMAP^ 2 = 0.1099 from WMAP data [104]. For the error we have chosen a nh 2 = 0.01 
to allow for some theoretical error on top of the experimental one. We have also chosen 
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importance functions where we add factors to drive the scans into regions of high gaugino 
fractions, high scattering cross sections or low neutralino masses. To get a well-sampled 
database we have also forced ADSCAN into certain mass ranges with importance functions 
of the type 



g (n x h 2 .. 



m x o ) = exp 



1 fn x h 2 - HwMAP h 2 \ 2 \ ( m X \~ target N 2 



2 V a Qh 2 J 2 V Am 



varying mtarget between 10 GeV and several TeV, and commonly using Am = 20 GeV. 

In Fig. 6 we show the results of a model exclusion analysis performed on this database 
for an 86-string IceCube configuration (IC86). 

The detector will be most sensitive to the region where the neutralino is an almost 
equal Higgsino-gaugino mixture, which is where the spin-dependent scattering cross-section is 
usually the largest (as can be seen especially in the right panel in Fig. 6). This occurs because 
the spin-dependent cross-section is dominated by t-channel Z exchange, which couples to the 
difference between the two Higgsino contributions; if the neutralino is strongly gaugino, this 
coupling is absent, whereas if it is strongly Higgsino it becomes equal parts Hf and H® an d 
the two contributions cancel. 

As many of these models are captured in the Sun mainly via spin-dependent scattering, 
the current direct detection experiments (which are mostly sensitive to the spin-independent 
scattering cross section) do not probe most of the models to which IC86 will be most sensitive. 
This can be seen in the middle panel, where many of the models IceCube can probe are far 
below the sensitivities of both current and near-future direct detection experiments. 

Compared to the standard IceCube 'hard' and 'soft'-channel analyses [41], which are 
valid for just two specific annihilation final states, the (mixed final state) MSSM-7 models 
excluded in the analysis here lie (predictably) between the hard and soft sensitivity curves. 
The implied sensitivity is closer to the soft-channel curve than it would have been if we 
had employed a similarly aggressive angular cut as in [41]. For this analysis we used an 
angular cut of (f> cu t = 20° around the solar position. In general we have chosen angular 
cuts so as to include the majority of a potential signal from the Sun. The cut angle may be 
optimised for different models however, as higher-mass WIMPs would certainly produce more 
centrally-localised, high-energy neutrino events. This is particularly important when doing 
model exclusion, as the p value Eq. 3.28 depends only on the expected number of signal and 
background events, not their angular distributions. The exclusion power of the full likelihood 
Eq. 3.20 is less sensitive to the cut angle, but is still not independent of it, because of the 
necessary approximations we had to make about the separability of the spectral and angular 
components. In principle the cut angle could be profiled or marginalised over as a nuisance 
parameter in model scans, but this is probably overly conservative; it can be optimised 
and chosen in advance for specific models based on simulated data, and an approximate 
scheme (depending, say, on WIMP mass and annihilation branching fractions) constructed 
for choosing the approximate optimal cut angle for any arbitrary model. Such an extension 
may be implemented in future updates to the formalism we have presented in this paper. 

5 Conclusions 

We have constructed an unbinned likelihood formalism for including full event-level IceCube 
data in parameter explorations of theories for new physics. The likelihood function includes 
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information about the number, direction and spectral characteristics of neutrino events, and 
is fast to calculate. We also constructed a simple associated measure for model exclusion, 
which is even faster to calculate and can be used for single models, without any reference to 
other parts of a parameter space. 

We performed a number of example model scans using our likelihood construction within 
the MSSM. We carried out global fits to the CMSSM with actual 22-string data, and with 
the 22-string effective area rescaled to represent the final detector configuration. Existing 
22-string data has little impact on the CMSSM, but the final detector configuration will 
robustly exclude the majority of the focus point region if IceCube finds no evidence for 
WIMP annihilation. We carried out a mock signal recovery in the CMSSM, showing that 
our method accurately recovers a benchmark point, and constrains model parameters very 
well. In the process, we showed the utility of including spectral information in IceCube 
searches for dark matter. Finally, we carried out a simple example model-exclusion analysis 
in the MSSM-7, showing that the 86-string configuration of IceCube will test a number of 
models that cannot be probed by direct detection experiments in the near future. 

Our likelihood construction is implemented in DarkSUSY 5.0.6, and freely accessible to 
the community. It will also be made available in a future release of SuperBayeS. The data 
and simulations we have used for this study, including 22-string IceCube event lists [24], are 
freely available on the web [80], and in DarkSUSY 5.0.6. 
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